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In the bistable regime of the FitzHugh-Nagumo model of reaction-diffusion systems, spatially homo- 
geneous patterns may be nonlinearly unstable to the formation of compact "localized states." The 
formation of space-filling patterns from instabilities of such structures is studied in the context of 
a nonlocal contour dynamics model for the evolution of boundaries between high and low concen- 
trations of the activator. An earlier heuristic derivation [D.M. Petrich and R.E. Goldstein, Phys. 
Rev. Lett. 72, 1120 (1994)] is made more systematic by an asymptotic analysis appropriate to 
the limits of fast inhibition, sharp activator interfaces and small asymmetry in the bistable minima. 
The resulting contour dynamics is temporally local, with the normal component of the velocity 
involving a local contribution linear in the interface curvature and a nonlocal component having the 
form of a screened Biot-Savart interaction. The amplitude of the nonlocal interaction is set by the 
activator-inhibitor coupling and controls the "lateral inhibition" responsible for the destabilization 
of localized structures such as spots and stripes, and the repulsion of nearby interfaces in the later 
stages of those instabilities. The phenomenology of pattern formation exhibited by the contour 
dynamics is consistent with that seen by Lee, McCormick, Ouyang, and Swinney in experiments on 
the iodide-ferrocyanide-sulfite reaction in a gel reactor. Extensive numerical studies of the underly- 
ing partial differential equations are presented and compared in detail with the contour dynamics. 
The similarity of these phenomena (and their mathematical description) with those observed in 
amphiphilic monolayers, Type-I superconductors in the intermediate state, and magnetic fluids in 
Hele-Shaw geometry are emphasized. 



I. INTRODUCTION 

Recent experimental studies of pattern forma- 

tion in reaction-diffusion systems have revealed a mech- 
anism for the generation of space-filling patterns that 
is markedly different from the classical Turing bifur- 
cation 1^. In the Turing scenario, realized recently 
in several experiments [^^|^, a periodic pattern arises 
throughout space from a linear instability of a homo- 
geneous state. In contrast, experiments of Lee, Mc- 
Cormick, Ouyang, and Swinney have shown that 
space-filling "labyrinthine" patterns can develop from 
finite-amplitude perturbations to linearly stable homo- 
geneous states (Fig. |^). These labyrinths are charac- 
terized by patches having different chemical composi- 
tion that are separated by relatively sharp interfaces, or 
fronts. The pattern formation process involves the mo- 
tion of these fronts, which are generally observed to be 
mutually repelling, thus preventing self-crossings and as- 
sociated changes in topology. The observations that this 
system possesses bistability, requires finite-amplitude dis- 
turbances for nucleation of patterns, and displays fin- 
gering instabilities of compact domains suggests a con- 
nection to the one-dimensional, reaction-diffusion "local- 
ized states" considered by Koga and Kuramoto . Sub- 



sequent generalizations to higher dimensions by Ohta, 
Mimura, and Kobayashi |^ revealed that these localized 
states could exhibit fingering instabilities |^ . 

In earlier work |lO) it was suggested on the basis of 
heuristic arguments that this kind of pattern forma- 
tion by interacting chemical fronts could be understood 
by means of a "nonlocal contour dynamics model," de- 
rived from the well-known FitzHugh-Nagumo model of 
activator-inhibitor competition |liyi2f| in the limit of fast 
inhibition. This law of motion derives from the combina- 
tion of a Young-Laplace force associated with interface 
curvature [ p^ and a screened Biot-Savart coupling be- 
tween distant segments of the interface. This nonlocal 
contribution embodies the phenomenon of "lateral inhi- 
bition" - an inhibitory action on a scale much larger than 
the activator front thickness [Q^. Numerical studies 
of the contour dynamics model revealed the essential fea- 
tures seen in experiments: the instabilities of compact 
structures, repulsion of chemical fronts, and the relax- 
ation of branched structures to compact ones as control 
parameters are varied. 

The competition between Young-Laplace and Biot- 
Savart forces has been shown to appear in a variety of 
other pattern-forming systems, each of which exhibits 
labyrinthine interface evolution. These include monolay- 
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ers of dipolar molecules at the air- water interface [ [l6[ , 
magnetic fluids in Hele-Shaw flow |l8| Jig|,|20| , Type-I 
superconductors in the intermediate state and thin 
garnet films ||2^. In each of these cases, the patterns 
are defined by the boundaries between two thermody- 
namic phases in coexistence, with bulk electric or mag- 
netic dipolar order. The Biot-Savart interactions then 
arise from the usual correspondence between magnetiza- 
tion and current loops. The common phenomenology of 
these systems has been reviewed recently 

Here we elaborate on the contour dynamics in two 
ways: a systematic derivation is presented using matched 
asymptotics, and extensive numerical evidence is pro- 
duced to confirm consistency with the phenomenology of 
the original reaction/diffusion partial differential equa- 
tions (PDEs). Throughout the analysis we emphasize 
the variational structure of the dynamics, both in the 
fast-inhibitor limit and more generally, for it appears not 
to be widely appreciated that complex patterns can arise 
from purely gradient flows. 

In order to motivate the major emphases in this work, 
we show in Figures |^ and ^ two numerical simulations 
of the pattern formation exhibited by the FitzHugh- 
Nagumo PDF model. The regions in which the activator 
takes on each of its bistable values are indicated by black 
and white. In Figure ^ we see a compact initial domain 
of one phase expand and finger to produce a space-filling 
labyrinthine pattern. Yet a small change in parameter 
values leads to the relaxation of this pattern back to a 
stable localized state, as shown in Figure |. In this paper 
we show that the complex pattern formation observed in 
Figures ^ and |^ can in fact be explained in terms of three 
elementary processes; 

• stripe stabilization and repulsion 

• domain localization 

• interface proliferation and transverse front instabili- 
ties 

The remainder of the paper is organized around the eluci- 
dation of each of these features on both energetic grounds 
and by means of asymptotic methods for the derivation 
of front motion from PDFs. The particular reaction- 
diffusion model of interest here is introduced in Section 
II along with a discussion of its stability and variational 
structure. We illustrate the region in parameter space 
where the Turing instability is precluded and localized 
states appear and discuss the gradient-flow nature of the 
dynamics in the fast-inhibitor limit. Interface dynamics 
in this limit is considered in detail in Sections III and 
IV, where we reiterate the heuristic arguments leading 
to front dynamics and present the asymptotic analysis. 
Section III focuses on one-dimensional systems where the 
features of stripe stabilization (a form of domain local- 
ization) and stripe repulsion are most easily understood. 
Two-dimensional systems are considered in Section IV, 
in which a variant on domain localization is found. A de- 
tailed discussion of transverse front instabilities is given 
in Sec. V, revealing a common mechanism not unlike that 
of the MuUins-Sekerka instability in solidification. The 



numerical studies of both the full PDEs and the contour 
dynamics are presented in Section VI. Some considera- 
tions regarding the form of the dynamics away from the 
fast-inhibitor limit are outlined in Section VII. Nonlo- 
cal interface dynamics in the various contexts described 
above are tied together in Section VIII, and Section IX 
outlines our conclusions and open problems. For com- 
pleteness, an appendix details the numerical methods 
used both for the study of the PDEs and the contour 
dynamics. 

There have been several other recent theoretical stud- 
ies relevant to the experiments of Lee, et al. In ex- 
tensive simulations of the Gray-Scott model |^^, Pear- 
son 1^ has found a wide variety of patterns, includ- 
ing ordered arrays of spots, lamellar stripe domains, and 
labyrinthine structures. Some of these were found to 
arise from finite-amplitude perturbations, as in the ex- 
perimental work. More recently, Meron and Hagberg 
P6| , p7| , based on reaction-diffusion equations similar to 
the FitzHugh-Nagumo model, have provided elegant geo- 
metrical arguments for a connection between labyrinthine 
instabilities in the fast-inhibitor limit and spiral wave be- 
havior for slow inhibition. The present work is comple- 
mentary to these in providing a more detailed picture of 
the fast-inhibitor limit, where the dynamics reduces to 
a pure gradient flow. Clearly, however, certain phenom- 
ena are precluded by the fast-inhibitor limit, such as spot 
"self-replication" studied in recent experimental |2^ ] and 
theoretical works j2^j2^,^ . 

II. THE MODEL 

A. Definition 

Our starting point is the FitzHugh-Nagumo model |Tl| ] 
for the coupled dynamics of an activator u and inhibitor 
V. We choose appropriate definitions of u, u, space and 
time to arrive at the nondimensionalized form 

ut= DV'^u- F'{u]r)- p{v-u) (2.1a) 

evt = V^w -v + u . (2.1b) 

Here, D is the activator diffusion constant normalized to 
that of the inhibitor, and F{u;r) is a double-well poten- 
tial representing the autocatalytic behavior of the acti- 
vator. In a convenient parametrization whereby the two 
local minima in F are at u = and u = 1, the potential 
F{u; r) is written in terms of a single symmetry-breaking 
control parameter r as 

F(u:r) = -u'^ (u- i f + ( r - -] ( -u'^ - -u^ - — | 

(2.2) 

provided < r < 1. When r differs from 1/2, a difference 
in potential is created between the two states. 
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AF = F{u = 1; r) - F{u = 0; r) 
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(2.3) 



so that for AF > (r > 1/2) the state m = is the 
more stable, and the reverse for AF < 0. In Figures ^ 
and ^. the u = state is associated with the white area 
and u = I with black. 

With the potential F as defined above, the model PDE 
is a coupled diffusion system which is nonlinear only in 
the activator equation 



Ut 



DV'^u- u{u- r) {u- I) - p{v - u) (2.4a) 



V^w — [v — u) 



(2.4b) 



This particular nondimensional formulation has the ad- 
vantage of possessing an exact invariance under the si- 
multaneous transformations 



l-u 



1 - V 



l-r . 



(2.5) 



This property implies that we need not consider as dis- 
tinct cases the evolution of black spots in white domains 
versus white spots in bl ack domains. Within the con- 
text of the PDE model (2^) we therefore focus on the 
question: 

"By what mechanism do black labyrinthine patterns 
develop in a white domain?" 

The form of the inhibitor dynamics embodies both its 
self-limiting behavior and its stimulation by the activa- 
tor, while the linear coupling to v in the activator dynam- 
ics is the simplest such term. Although it is traditional 
to view V itself as the inhibitory agent, we have chosen 
instead to introduce an inhibitor coupling in the form 
p{v — u) which insures that the state u — v — 1 remains 
a stationary uniform state fo r all p. 

The parameter e in Eq. 2.1b distinguishes between 
the slow-inhibitor limit (e ^ 1) and the fast-inhibitor 
regime (e ^ 1). We assume the latter, a limit that is 
opposite to the limit assumed in phase-field models ||3l| ] 
and spiral wave dynamics | |3^ . Through the nondimen- 
sionalizations used in (2.4), the natural diffusive length 
scale of the inhibitor is set to unity. The corresponding 
length scale for the activator is a/D, so the limit of sharp 
activator fronts of interest here requires ^/D <C 1. The 
parameter D serves to mainly to set the width of those 
fronts, while the two remaining parameters in the model, 
r and p, will be seen to control more global aspects of 
the front dynamics. 



B. Stability against the Turing bifurcation 

A simple stability analysis shows that both states u = 
and u = 1 can be linearly stable simultaneously to 



all periodic disturbances. Consider first the state u = 
V = 0. In the limit of fast inhibition, e — > and for 
< r < 1 (so that both minima of F exist), one finds 
that perturbations having the form 



(2.6) 



are characterized by two branches of solutions for the 
growth rates cr. 



CT+ = -Dk^ - ■ 
1 + P 



fc2 



' H 
0{l) 



0{e) 



(2.7) 



While the branch ct- is clearly damped for all wavevec- 
tors, (7+ may become positive for p exceeding a critical 
value Pc{t)- The neutral curve for this linear (Turing) 
instability, as defined by the conditions <J+{k) =^ and 
da+{k)/dk ~ 0, is 



prir) ^ (^V^+Voy 



(2., 




so that when p < PTi^) the state u = u = is linearly 
stable. At the critical p = prir) the marginally stable 
wavevector kx is 



(2.9) 



Appealing to the symmetry relations in Eq. p.5| , we 
deduce that the stability criterion for the state u = v = 1 
i s p < pt(1 — ''), where the function pT the same as in 
( ^.8| ). Thus, simultaneous stability of both homogeneous 
states to the Turing bifurcation requires p < pt{t) and 
p < pt{^ — r). The region in p — r space so defined 
is shown in Figure || for the case D = 0.01. The inset 
to Fig. ^ shows two curves of a+{k) for the parameters 
(r = 0.65, p = 0.25) (stable) and (r = 0.65, p = 0.60) 
(unstable). 

As the markers on Figure ^ suggest, the labyrinthine 
dynamics of Figures I and I are observed near p « 
and r « 1/2 - values which are well within the sta- 
ble regime. Although for these parameters the system is 
stable against the Turing mechanism, it is the underly- 
ing bistability that permits the existence of non-trivial 
localized states through intrinsically nonlinear processes. 



C. Variational Aspects 



Equations 2A are neither purely dissipative nor Hamil- 
tonian, but can be expressed in a variational form 



Ut = 



5£u bT_ 

5u ^ 5u 



(2.10a) 
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Sv Sv 

where the energy functionals are 



(2.10b) 



Jd^l^^D\Vu\'' + F{u;r)-^pu^^ 



T= d: 



(2.11) 



It is precisely because the cross-terms in the dynamics 
(2.4) are of opposite sign that the system is not a pure 
gradient flow for finite e. One may verify that neither 
nor f nor T decreases monotonicafly in time. This 
completely chan ges in the fast-inhibitor limit (e 



0), for then Eq. 



2.10b becomes a functional relation 



5£v 8T 
5v Sv 



(2.12) 



from which it follows that the combination + pj- de- 
creases monotonically in time 



d_ 

di 



{£u + pJ') = -Jd' 



5£ ST 
X I —+p — 

ou Su 



(2.13) 



The dynamics is then a gradient flow, a feature central 
to the heuristic arguments given in earlier work and 
elaborated on below. 



D. The Nonlocal Energy Functional 

With the exception of Section VII, for the remainder 
of this paper we shall assume the fast-inhibitor limit 
e = 0. The inhibitor dynamics degenerates into an 
inst antan eous-in-time relation between v and u given by 
Eq. 2.12 . Using the specific forms of £v and Ti it reduces 
to 



(V^ - 1) t;(x,t) = -u(x,t) 
which can be solved using a Green's function 

v{'K.,t) — j dx'Q{x ~ x')u{x') . 



(2.14) 



(2.15) 



In one and two dimensions respectively, the Green's func- 
tions are 



G{x — x') — \e 



1 



(2.16a) 



Cy(x-x') = — i^o(|x-x'| 



(2.16b) 

where i^o is the modified Bessel function of order zero. 



Substituting for v in Eq. 2.4a. we obtain the spatially 
nonlocal activator PDE 



Ut 



D\/'^u - F'{u; r)+ pu- p J dx'g{x - x')u(x') . 



(2.17) 



From the arguments leading up to Eq. 2.13 , and as had 
been observed earlier by Ohta, Ito, and Tetsuka [Q, Eq. 
2.17 has the variational form ut = ~5£ /Su with the en- 
ergy functional 

+ ^p JdxJ dx'M(x)g(x - x')m(x') . (2.18) 

While, as remarked earlier, £ decreases monotonically 
in time, £ is not necessarily bounded from below in an 
infinite domain, and this variational principle does not 
guarantee a stationary long-time limit. 

As an aside, if the variations in u are on length scales 
long compared to the 0(1) inhibitor screening length, 
then we can expand the integrand in ( ^.18 ) for x' near x 
and obtain up to fourth order the gradient expansion 



(2.19) 



Collecting together terms, we obtain an approximate lo- 
cal theory, 

£c^ Jdxi^^D\Vu\' + ^p{\7\f +Fiu;r)^ (2.20) 

where D = D — p. At this level, the energy func- 
tional is very similar to that which appears in the Swift- 
Hohenberg model as well as the theory of Lifshitz 
points in condensed matter systems , and leads quite 
naturally to modulated patterns by virtue of the possibly 
negative coefficient of | Vup. We note however, that such 
a gradient expansion is only of limited utility in the limit 
of sharp interfaces of interest in the present work. Nev- 
ertheless, an approximate local theory can be developed 
within the limit of sharp interfaces when the curvature 
of the interface is small, as shown in Section V.B. 



Let us note finally one important feature of Eq. 2.2C. 
We s ee tha t the lowest-order term in the gradient expan- 
sion (2.19) leads to a contribution {l/2)pu^ in the energy 
functional. This precisely cancels the term — (1/2) pu^ 
seen in Eq. 2.18 which arises from the —p(v — u) inhi- 
bition coupling. Thus, the effective potential in a local 
theory (2.20) is precisely the function F{u;r), a conve- 
nient simplification in subsequent sections. 
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III. INTERFACE MOTION IN ONE DIMENSION 
A. Energetic Derivation of a Stripe Evolution 

In the PDE dynamics shown in Figs. ^ and ||, the po- 
tential difference AF is positive so that u = (white) 
is the energetically-preferred state. Nonetheless, in both 
cases the black domain survives despite its energetic dis- 
advantage. Insight into this bistable configuration can 
be obtained through an analysis of a single black stripe. 

Consider the dynamics of an interval with u w 1 in a 
background of u w (i.e. the c ross -section of a black 
stripe). A PDE computation of ( ^.4|) for one such evo- 
lution is shown in Fig. O, where we see an approach of 
the two sharp activator Ju) fronts that mark the transi- 
tion from black to white. In contrast, the inhibitor field 
V is considerably more diffuse. But, perhaps the most 
notable aspect about this particular evolution is that the 
black stripe asymptotically stabilizes at a finite size. By 
applying a simple heuristic argument that was outlined 
in an earlier work |lO| , an approximate dynamics for this 
process can be derived which predicts the equilibrium 
width of this black "localized state." 

The argument begins from the essential assumption 
that the activator fronts are narrow relative to the scale 
of the pattern, so that u may be taken to be piecewise 
constant away from fronts. We focus on the case of a 
single black stripe of u = 1 located symmetrically be- 
tween X — —Qit) and X — +Q{t) in an otherwise quies- 
cent white background {u = 0). The evolution of Q{t) is 
deduced from the gradient dynamics based o n th e one- 
dimensional version of the energy functional ( p. 18 ) 



£[u] 



1 2 

dx < —Du^ 



F{u-r) 



-pw 



"2^ 



dx 



dx'u{x)g{x-x')u{x') . (3.1) 



To insure a finite energy, we consider the energy differ- 
ence Af = £[u\ —£[uoc\, where Uqo is the value u takes as 
X — *■ ±oo. Then, in following the prescription of the ear- 
lier results |To[ , the energy functional for a stripe is par- 
titioned into three contributions: an effective line tension 
7, a pressure 11, and a nonlocal coupling between fronts. 
These contributions are explicitly calculated below. 

In calculating A£ for the stripe, the line tension con- 
sists of the contributions that are localized at the fronts 
(like M^). In one dimension, this requires a contribution 
to A£ that is some constant 7 proportional to the num- 
ber of fronts. On dimensional grounds, the line tension 
associated with a single front 7 can be estimated by the 
gradient integral 



7 ; 



dxDui^O 



front 



D 



(3.2) 



where £ is a characteristic length scale of the front. Since 
this length scale must be the diffusive length a/D, we 
expect 7 ^ ^/15. 



Under the assumption that u is piecewise constant and 
either zero or unity, the evaluation of the pressure and 
nonlocal contributions to the energy functional are con- 
siderably simplified. By virtue of having subtracted the 
white background energy £[uoo\, the contribution arising 
from the potential F is restricted to the stripe, and is 
thus proportional to the pulse length 2Q. The constant 
of proportionality associated with the integral of the po- 
tential function F is interpreted as a pressure 11 and is 
given by 



^£'^dx{F{l;r)~F{0;r)}^AF. (3.3) 



The remaining two terms in the energy integral (3.1) 
are proportional to the inhibition coupling p and are also 
integrations only over domains where u — I. For the 
single black stripe, these two integrals, when combined, 
require the evaluation of 



dx<l- 



dx'Q[x — x') 



(3.4) 



Using the PDE for the Green's function Q ~ Q^x + 5{x — 
x') a nd the reciprocity relation Q^x — ~Qxx' , the integral 
(3.4) simplifies to 



1 

-p / dx / dx'Gxx' = p [Q{2Q) - g{Q)] 

^ J-Q J-Q 



(3.5) 



Combining all three contributions (3.2, p.3| , 3.5) the 
energy of a localized state of size 2Q is thus estimated to 
be 



A£[Q] ~ 27 + 2g n + p [g{2Q) - g{o)] 



(3.6) 



As remarked in the paragraph following Eq. 2.2C , we see 



that the pressure contribution 11 comes from the bare 
function F(u;r), rather than from the shifted function 
F- (l/2)pit2 in Eq. 



2.18. Furthermore, the line tension 



constant 7 as estimated from Eq. 3.2 is in fact asymptoti- 
cally correct, despite an apparent discrepancy of a factor 
of two between (3.1) and (3.2). This detail is resolved 
by the asymptotic analysis in the following section, and 
involves the recognition of an additional 0{^/D) contri- 
bution from a boundary correction to the estimate in 

O- 

The final step in this heuristic derivation is the trans- 
lation of the variational principle ut = ~5£ /5u into an 
equation of motion for the front position Q{t). In terms 
of the general theory of Lagrangian dynamics, viscous 
forces can be introduced as the variational of a Rayleigh 
dissipation functional TZ p7|] , 



1 /■+°° 
^ J -00 



(3.7) 



so that the dynamics can be written in the variational 
form 



5 



d dC 9£ 

dt dut du 



dn 

dut 



(3.8) 



For a single moving front having an invariant spatial pro- 
file, the dissipation functional is intimately related to the 
line tension constant 7 since 



n[ut] 



dx 



front 



D 



(3.9) 



where the result is doubled to account for both fronts. 
Now, adopting Q as the dynamical coordinate and taking 
for the Lagrangian the en ergy difference C[Q] — — A£[(3], 
the variational principle (S.8|) gives a dynamical equation 
for the front location Q{t) 



Qt = - 



D OAS 



D 

7 



2 



(3.10) 



Note that the absence of any ut-dependence in the La- 
grangian £[u] leads to an inertialess front dynamics. The 
possibility of inertial contributions to the front evolution 
is discussed Section VII. 



We deduce from ( 3.10 ) that black localized states are 
only possible if AF > 0, so that the state u = 1 has 
the higher energy of the two minima of i^(M). Moreover, 
beyond the threshold p — 2AF the inhibitor coupling 
supports a stable black stripe of width 



2Q* ^ In 

and whose associated energy is 
2A£ V2 



r 
3^ 



1 - In 



3p 



(3.11) 



(3.12) 



We have introduced the scaled parameters 
, _ 6AF _ r - 1/2 , p 



D 



D 



P=^, (3.13) 



which in the following section are shown to follow nat- 
urally from an asymptotic analysis. Below threshold, 
the inhibition coupling is insufficient to prevent annihi- 
lation through the collision of two fronts. On the other 
hand, for AF < 0, no equilibrium width exists and the 
stripe expands indefinitely to fill space. The existence 
of an equilibrium stripe width for the less-preferred state 
demonstrates, within the context of this PDE model, the 
process of stripe stabilization. 



Using the symmetry relations (2_^), we deduce the 
converse result that white localized states {u = 0) in a 
black background {u = 1) may exist only when AF < 0. 
Likewise, a white stripe will expand unboundedly when 
AF > 0. 

Localized states may also be formed by the nucleation 
of a small domain of u = 1, followed by its growth to 



the stable size. From the form of the energy (3.6) we see 
that there is a barrier to the creation of a localized state, 
so this process requires a nucleus larger than some criti- 
cal size. This is illustrated by the PDE computation in 
Fig. 1^, where we see the growth of a slightly supercritical 
nucleus. 

While we have emphasized the somewhat counterintu- 
itive result that a localized one- dimensional state consists 
of the phase which has a high er va lue of the potential F, 
it is interesting to note from ( |3.12| ) that the total "struc- 
tural energy" A£ of the stripe can be of either sign. It 
will ultimately be shown (Sec. ^ that the existence of 
stripes with AS = is intimately connected with the 
onset of the labyrinthine instability. 



B. Asymptotic Derivation of a Front Dynamics 

The previous section describes the dynamics and inter- 
action of interfacial fronts as a Lagrangian relaxation pro- 
cess which is independent of the detailed fine structure of 
the activato r (u) fi eld. Using formal matched-asymptotic 
arguments ||l^,^,^,|8| , these results are substantiated 



as the leading order behavior of the PDE model 2.1 in a 
particular limit of weak activator diffusion. 

The asymptotic analysis relies on a strong separation 
of scales between the spatial structures of the activator 
and inhibitor fields. Through t he n on-dimensionalization 
of the reaction/diffusion PDE 2.1, the 0(1) spatial scale 
is defined to be the natural diffusion length of the in- 
hibitor (v) field. In the limit of weak activator diffu- 
sion (D <C 1) the transitions of the activator field be- 
tween bistable states are resolved within thin interfaces, 
or interior layers, whose 0{Vd) width is characteris- 
tic of diffusion-limited fronts. From the perspective of 
matched-asymptotics, the sharp activator fronts on the 
inner scale are contrasted by the more gradual relaxation 
of the inhibitor concentration that occurs on the outer 
scale. 

Outer 0(1) Scale: In one dimension (x), the model 
equations 2.1 are rewritten explicitly 



u{u - r){u - 1) = Duxx - piv - u) -ut 



{v- 



evt . 



(3.14a) 



(3.14b) 



where the activator diffusion {D <^ 1) and inhibition 
coupling (p <C 1) are assumed small. In anticipation of 
a dynamics consisting of slow front motion (implied by 
D 1), both time-derivatives may also be co nsistently 
assumed small. The right-hand sides of 3.14 are then 
seen to be perturbations with respect to this outer scale. 

In terms of the outer variables, the underlying pattern 
is defined by regions where the activator is in eith er of it s 
locally stable values. For the cubic nonlinearity (3.14a), 
these are 
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which defines a binary pattern on the line. The present 
approach is more general than the heuristic analysis of 
the previous section, capable of treating a pattern con- 
sisting of any number of fronts. The leading order outer 
inhibition field, v^{x), then satisfies the inhomogeneous 
elliptic equation 



1 black areas 
white areas 



(3.16) 



where continuity of v'^ and is imposed at jump discon- 
tinuities in the binary pattern u'^{x). As well, suitable 
boundary conditions at the domain edges (for example: 
periodic, or = at infinity) must be included to insure 
a unique solution. 

The time dependence of this leading order outer so- 
lution is implicit in the functions Q-'{t) that locate the 
boundary points between the black and white areas, the 
jump transitions of vP . The velocity Ql of these inter- 
face points is obtained as a consequence of the asymptotic 
matching at inner scales to remove the discontinuities of 
the chemical fields. 

Inner 0{\/lj) Scale: At each interface point x — Q^{t) 
the discontinuity of the outer activator solution vP{x) is 
resolved by an inner representation defined locally on a 
finer diffusion scale rj = {x — Q^)/\/D. Since the activa- 
tor solution is not discontinuous on this inner scale, the 
precise location of the interface {rj — 0) is defined by the 
condition u{x = Q^) = 1/2. In terms of this rescaled 
and moving coordinate, the equations for the inner so- 
lutions U{ri,t) and V{i],t) reflect a markedly different 
asymptotic ordering. Introducing the nonlinear operator 
S[U] 



S[U]=U„„^U {U- 1/2) (U-l), 
the inner equations become 

S[U] = --^Un - 6AFU (U - 1) 



(3.17) 



+p iV-U) + Ut 



K,r, ^D{V-U) + eDVt 



(3.18a) 



(3.18b) 



The homogeneous solution of the left-hand sides of ( 3.18| ) 
correspond to a stationary front solution 



1 ^ tanh 



\2y/2 



V° = v^{x = Q^) , 



(3.19) 



which, as 77 ^ ±oo, match to the behaviors of the outer 
solutions M°(a;) and v^{x) at the interface x — > {Q-')^- 
The sign =F is chosen by the orientation of the front, 
where a negative sign corresponds to a left-to-right tran- 
sition from w = 1 to u = 0. Note also that, to leading 
order, the inhibitor field is constant within this interior 



layer - this is co nsiste nt with the continuity of the outer 
solution w°(x) of ( 3.16| ) which only requires an 0{D) tran- 
sition solution to resolve the discontinuity in its second 
derivative. 



The stationarity of this leading order front ( 3.19 ) is 
broken by the effects of the right-side terms of ( 3. 181) and 



suggests a perturbative derivation for the front speed Q/. 
based on the asymptotic balance 



Ql 



AF - p < 1 



(3.20) 



The determination of Q/, w hich is essentially just a non- 



linear eigenvalue for ( |3.18 ), then follows from the solv- 
ability condition for a (bounded) first correction to the 
inner activator solution 



C/(r?)~f/°(??) + pf/^(??) . 



(3.21) 



Substitution into ( [3.18a| ) gives an inhomogeneous ODE 
for the time-independent first correction U^^ij) involving 
the operator 5'[f/"'^] tha t is l ust the linearization of the 
nonlocal operator S[U] ( 3.17 ) 



^ rjri 



F"(C/";l/2) C/i 



(3.22) 



where the primes on F denote differentiation on the first 
variable. The first-order perturb ation s are restricted to 
the terms involving the balance ( [3.20| ) and gives 



6AF 



2D 

+p [v\x^Q')-U'' 



C/° (t/" - 1) 



(3.23) 



Application of the identity ^2 [7° = =fJ7° (f/°- 1) has re- 
sulted in some simplification of terms in the above equa- 
tion. 

By the translation symmetry associated with the 
[/"-solution, the app ropri ate solvability condition for 
bounded solutions to (3.23) is the inner product integral 



(3.24) 



Since the linearized operator S' is o rthog onal to 
the inner product of the equation ( 3.23 ) determines 
the unique front speed that ensures the existence of a 
bounded first correction U^{r]) 




(3.25) 



where again the choice of sign =p is consistent with the 
choice made for the leading-order solution (3.19). Self- 



consistently, the front velocity satisfies <C OVID. 



The important feature of the front speed result ( 3.25| ) 
is its dependence on the strength of the outer inhibition 
v'^{x = Q-')- By the elliptic nature of the outer inhibitor 
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equation ( ^.16 ), v{x = Q^) will in general be influenced 
by the location of other fronts - and hence embody the 
nonlocalitV; Thus for a sequence of fronts, the speed 



relations ( 3.25 ) for j = 1 to will be a set of N coupled 
ODEs that describe the boundaries of the binary pattern 
for u{x). 

Motion of a single front. For this case, on the infinite 
line, we consider a pattern that is white (u = 0) to the 
right of an interface at a; = Q{t) and black (w 1) to the 
left, 



u°{x) 



1 x< Q{t) 
Q{t) < X 



(3.26) 



Solving (3.16), the outer inhibitor field has the simple 
piecewise continuous solution 



v°(x) 



X < Q{t) 
Q{t) < X , 



(3.27) 



where — boundary conditions are imposed at infin- 
ity. Since v''^{x = Q) = 1/2 at the interface, by ( 3.25 ) 
there is seen to be no 0{p) correction to the front speed 



Qt 



(3.28) 



so that, when AF > 0, the front travels with constant 
negative (leftward) velocity. Motion towards the black 
region is the expected result when the white [u — 0) 
state is energetically-preferred (^). 

Higher-order corrections to the speed of a single front. 
The calculation of the front profile and speed may be 
continued to another order in p in the case of a solitary 
front. The solution of ( 3.23| ) for the first correction is 



d 
drj 



Tytanh 



2V2 



(3.29) 



This describes the overshoot of the profile seen in Figs. 
^ and ^. The gradient of the front is increased in the 
neighborhood of 77 = 0, and this would be expected to 
increase the rate of dissipation, and hence retard the in- 
terface motion. This intuition is confirmed asymptoti- 
cally when AF ^ p ^ ^fO, as the next order of the 
solvability condition can be explicitly evaluated, giving 
the extended velocity condition 



[1 - 6p] 



(3.30) 



This linear correction to the front speed is clearly appar- 
ent in the plot of velocity versus p in Figure obtained 
by direct numerical solution of the PDEs. Note also that 
beyond p ~ 0.02, even this next order asymptotic correc- 
tion is insufficient for quantitative accuracy. 

Interaction of two fronts. Reconsider the situation in 
an infinite domain —00 < x < +00 in which two inter- 
faces are symmetrically located at x = ±(5(i) between 
which u = 0. These fronts are the edges of the local- 
ized state considered in the heuristic derivation which 



obtained the dynamics ( 3.10| ) - this result is confirmed 
below using the matched-asymptotic front speed formula 

(iH). 

The calculation proceeds as in the case of one front, 
but begins from the outer activator pattern 



rO -00 <x<~Q{t) 
u"{x) = I 1 -Q{t) <x<+Q{t) 
lo +Q{t)<x <+oo 



(3.31) 



that is symmetric about the origin. Solving ( p^ , the 
outer inhibitor field can again be expressed in terms of 
exponential functions 



w°(a;) = 



sinh Q e^^ 

1 — e^'5 coshx 

sinh Q e~^ 



- 00 < X < ~Q{t) 

- Q{t) <x< +Q(t) 

+ Q{t) <x< +00 . 



(3.32) 



This determines the inhibitor contribution to the front 
speed 



.°(x = Q) = i(l-e-2Q) 



(3.33) 



where now, unlike in the case of a single front, v'^{x — 
Q^) 7^ 1/2 at the interface. Substitution into the asymp- 
totic formula for the front velocity (3.25) gives 



-2Q 



(3.34) 



where the min us sig n must be taken since u{x = Q~) — 1. 

The result (3.34) not only recovers the heuristically- 
derived stripe dynamics ( 3.1C| ), but also determines the 
constan t of proportionality associated with the line ten- 
sion 7 ( |3.2[) and the Rayleigh dissipati on ra te ( 3.9). Di - 
rect comparison of the two formulae (3.10) and (3.34) 
gives 
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(3.35) 



which is, as expected, 0(\/D). Moreover, since the line 
tension coefficient (^^) represents the energetic contri- 
bution proportional to perimeter, the leading order form 
of the front profile can also be used in the integration 
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(3.36) 



Note that the additional contribution of F{U ; 1/2) is a 
boundary error associated with the pressure integral ( |3.3| ) 
and resolves the discrepancy involvin g th e line tension 
energy that was alluded to below Eq. p^ . 

In addition to recovering the dynamics of stripe stabi- 
lization, the single front speed formula can be used to un- 
derstand the behavior of two black stripes when AF > 0. 
If we were to start from a configuration in which the 
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stripes were not at their equilibrium width, there would 
first be the ballistic collapse to form the black stripes of 
width near 2Q* . On a slower time scale, the two stripes 
repel each other by their exponential tails of inhibition 
(on the 0(1) scale) - this is an illustration of the process 
of stripe repulsion jTBI. 



C. Numerical Comparison 

Tests of the results for the dynamics and equilibrium 
separation of approaching fronts are shown in Fig. ^ for 
a situation like that in Fig. ^, but with the values D — 
0.0001 and p < 0.02, in the truly asymptotic regime. The 
front separation 2Q is shown as a funct ion o f time and is 
compared with the predictions of Eq. (3.34), using both 
the bare prefactor 6V2D associated with the leading- 
order profile, and the corrected prefactor 6^/2D{l~6p) of 
Eq. 
resu 



3.30. We see that the latter leads to highly accurate 



ts. The equilibrium separation between the fronts is 
independent of this prefactor, and as shown in Fig. ^ 
is in accurate agreement with the results of numerical 
solution of the full partial differential equations. 



IV. INTERFACE MOTION IN TWO 
DIMENSIONS 

A. Heuristic Derivation 

The heuristic derivation of the energetics and dynamics 
for contours in two dimensions proceeds along the lines 
used for one-dimensional patterns. This requires the ad- 
ditional assumption, beyond that of localized gradients, 
that the curvature of the boundary does not significantly 
affect the interface profile. We consider a single black do- 
main B of u — 1 surrounded by a sea of white, bounded 
by the contour r(s), where s is arclength. As in the 
calculation for one dimension, the energy functional AS 
consists of three contributions. 



The line tension contribution to the energy ( 2.1^ ) must 
now involve an integral around the boundary of the black 
domain 



dsj = 7L 



(4.1) 



dB 



where L is the perimeter, and where the line tension 
7 must be exactly that found earlier from the one- 
dimensional asymptotics (Eq. 3.35| ). 



By analogy with the result in Eq. 3.2, the potential 
integral contributes an energy proportional to the domain 
area A, 



I dx [F(l; r) - F(0, r)] = AFA . 
Jb 



(4.2) 



The main difference from the one-dimensional prob- 
lem is the treatment of the nonlocal term, which is pro- 
portional to the double integral / over the domain. 
Substituting Q — V^Q -f (5(x — x') we may recast the 
bulk contributions as boundary terms by twice using the 
Green's identity f . Vi/) = L ^ n?/; in a manner analogous 



to that in Eqs. 3.4 and 3.5 as follows 



dx I d-x! g - A = - i rfxVx • / dx'Vx'^ 
B JB Jb jb 

= - rfxVx • f ds'h{s')g 

JB JOB 

= - f ds (f ds'n(s) • h{s')g . (4.3) 

JdB JOB 

Note that n(s) is the unit normal vector pointing out of 
the black domain. Finally, collecting together all of the 
terms and noting that n(s) • n(s') = t(s) • t(s'), where t 
is the unit tangent vector, we obtain 

A£[r] =jL + AFA 

-ip <j)ds <j)ds'i{s) ■ i{s')g (r - r') . (4.4) 

Unlike in one dimension, the nonlocality cannot be re- 
duced to a pointwise evaluation of the Green's function. 

The form of the nonlocal coupling between tangent 
vectors is reminiscent of the self-induction of current- 
carrying wires, the direction of the current being spec- 
ified by t. Indeed, the sign of the interaction is such 
that antiparallel tangent vectors repel (like antiparallel 
current-carrying wires), while parallel ones attract. An 
alternative view of this connection with electr oma gnetic 
systems is to consider the inhibitor relation ( ^.14| ) as a 
Poisson equation relating a potential (w) to a charge den- 
sity (u), in which case the energy function J- — J uv in 
Eq. 2.11 is the associated electrostatic energy. Through 
the correspondence between magnetization and current 
loops, and the similarity of electric and magnetic dipolar 
phenomena, such a nonlocal interaction appears in the 
description of pattern formation in a variety of other sys- 
tems, as detailed in Section XIII. In such dipolar systems 
an important role is played by molecular or other cutoffs 
in the self-induction integrals. The present formulation 
neither has such cutoffs nor needs one for a well-defined 
energy, as the singularity in the Green's function g in 
(4^) as s — s' — > is integrable, being only logarithmic 
in |s — s'|. 

With the energy now formulated as a functional of 
the boundary contour r, it only remains to calculate the 
Rayleigh dissipation functional. This is just the integral 
around the boundary of the one-dimensional result 



n 
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2D 



f ds{n- TtY 

JdB 



(4.5) 



To obtain the functional derivative of the energy ([1.4|), 
we combine the local contributions 
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1 

— — = Kn 

V5 



(4.6) 



where g is the metric, with the result for a self-induction 
integral p8| , ^ 



j,ds (j)ds'i ■ i'g{R) 

ds'R(s, s') X i'g'{R) \ n(s) , (4.7) 



S_l 
6r2 



where R(s, s') = (r(s) — r(s')) /|r(s) — r(s')|, the cross 
product is a scalar in two dimensions (a x b = eijUibj), 
and the prime on Q indicates differentiation: Q' = 
dQ/dR — —Ki{R). The variational principle (3.8) de- 
termines the contour evolution, expressed as the normal 
velocity 



D 



n • rt = <^ 7'«(s) + 



p (pds'R{s,s') X i'g'{R) 



(4.8) 



Note that the nonlocal contribution in (4.8) is not a sin- 
gular integral, since the integrand for s' — > s has the 
finite limiting value 



jhn {R(s,5')xifi(|r(s)-r(s')|)} 



. (4.9) 



It is readily verified from Eq. 4.8 and Figure ^ that 
the vector product in the normal velocity embodies the 
repulsion between nearby antiparallel sections of the con- 
tour. As noted in earlier work [|l0|, this repulsion between 
adjacent fronts requires only that Q > and Q' < (for 
p > 0), and not on the details of Q. 



B. Asymptotic Derivation of a Curvature Dynamics 



Now we rederive the contour dynamics (4.8) by adapt- 
ing the asymptotic argument presented earlier for one 
spatial dimension. The only major modification involves 
the reorientation of the the inner coordinate (77) of the 
activator front to be in the normal direction to the inter- 
face. 

Outer 0(1) Scale: In two dimensions, x = {x,y), the 
leading-order outer representation satisfies the left-hand- 
side of the equations 

u{u - r){u - 1) = DV^u - p{v -u) -ut , (4.10a) 



V^v ~{v-u) = evt . (4.10b) 



The parameter scalings are chosen to obtain the analo- 
gous front dynamics, but with the added dimensionality, 
the activator diffusion scale (VD) is incorporated into 
the significant limit 



AF. 



y/D < 1 



(4.11) 



As before, the underlying pattern is defined by patches 
where the activator is in either of its bistable states 



u(x) 



black areas 
white areas . 



(4.12) 



Unlike in one-dimension where the transitions between 
black and white are points, in two-dimensions the inter- 
face separating a black patch from a white background is 
defined by a (slowly- moving) contour in the {x, ?/)-plane 



r(a;t) = {X{a:t),Y{a;t)) 



(4.13) 



where a is defined as a counter-clockwise (but not nec- 
essarily arclength) parameterization around the black 
patch. There, of course, may be more than one disjoint 
black patch in which case multiple contours r-' must be 
evolved - for simplicity of presentation however, these 
indices are omitted. The leading order inhibition field, 
v^{'x), then satisfies the inhomogeneous elliptic equation 
with a two-dimensional Laplacian 



black areas 
white areas 



(4.14) 



Along the interface r, continuity is imposed on u*^ and 
on (n • V)f'^, the normal derivative. Again, with suit- 
able boundary conditions at the domain edges, a unique 
solution is guaranteed. For consistency with the separa- 
tion of scales assumed in the asymptotics, it is essential 
that the interface geometry contain structure only on the 
inhibitor diffusion length - this demands that both the 
(global) interfacial separation distances and the (local) 
radius of curvature be 0(1). 

Local Orthogonal Coordinates The discontinuity of the 
activator on the outer scale is naturally resolved on an 
inner scale which is a stretched coordinate normal to the 
interface. At a point along the contour r{a;t), the local 
(unnormalized) tangent (t) and normal (n) directions are 
given by 



n = Jr„ 



+1 
-1 



(4.15) 
(4.16) 



where the matrix J corresponds to a 90° clockwise rota- 
tion. With the counter-clockwise orientation of the pa- 
rameterization a, the normal n points into regions of 
white. Using these as basis vectors for a local coordi- 
nate system, introduce a new coordinate fj that extends 
orthogonally from the contour r in the direction of the 
normal (n). This suggests the change of variable (sup- 
pressing the time-dependence) 
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x(ry, a) = r{a) + T]3ra{a) 



(4.17) 



so that fj ^ Q coincides with the contour r. In the thin 
interface hmit, the inner solution requires the vahdity of 
local coordinates only for small fj — 0{^/D). It proves 
calculationally more convenient to use instead an infinite 
(Lie) series as the change of variable 



x(77, a) = r + T]3Ya + ^'7^J^r, 
= exp SV-^-Q^ >^ 



(4.18) 
(4.19) 



for which (4.17) is the linearization. Straightforward dif- 
ferentiation demonstrates that the complete series (4.1E) 
has the property that it is orthogonal even for ?? 7^ off 
the generating contour r 



X|^ Jxq . 



(4.20) 



These are just the Cauchy-Riemann conditions, thereby 
showing that the change of variables between (x, y) and 
(ry, a) is locally conformal. Although for small fj only 
a few terms in this series are necessary in the present 
analysis, it establishes to all orders that the Laplacian 
nature of the diffusions are preserved up to a Jacobian 
metric. 



At each point along the interface r(a), the normal com- 
ponent of the velocity is sufficient to describe the inter- 
face motion. This freedom essentially follows from an 
arbitrariness in the contour parameterization a, which 
permits arbitrary choices of tangential velocity. Here we 
may choose the parameterization that yields, at fixed a, 
an interface moving normal to itself, with Qt{a;t) the 
leading-order normal velocity, 



(4.26) 



At the level of the first correction, the inner equations 
for th e fro nt profile are identical with the one-dimensional 
case ( 3.18 ) , but with the addition of curvature term from 
the Laplacian expansion (4.25) 



S[U] - --^Urj + '^VDktjU^,^ 



+p{V -U)+0{D) 



OiD) 



6AFU{U-l) 

(4.27a) 

(4.27b) 



where S is as defined above ( 3.17 ). The leading-order 
solution U^jri) is identical with the hyperbolic tangent 
solution ( 3.19 ). Application of the same solvability inner 
product (3.24) yields the front speed relation 



1 



|x„(r?,a;t)p 



(4.21) 



Inner 0{\'D) Scale: The resolution of the activator 
front, as well as the front speed determination now fol- 
lows analogously to the case of one dimension. The 
stretched inner coordinate rj 



1 



rj 



(4.22) 



resolves the activator front across the discontinuity, while 
the g-coordinat e lab els position along the contour. Using 
the expansion ( 4.19| ), we find that the metric factor in 
(Oil) is 



,(77,Q;;t)|" 



1 + 2y/DriK 



(4.23) 



where we have used the general expression for the curva- 
ture K, 



K = t • n. 



(4.24) 



defined such that it is positive for a convex black domain. 
Note that although the |rQ,|-normalization weakly breaks 
the orthogonality property, up to the first correction in 
curvature the Laplacian has the simple form 



V2 



hi 



1 - 2VDTJK 



92 



+ 0(1) 



(4.25) 



+P 



v°(a; = Q) 



(4.28) 



which clearly recovers the curvature and potenti al e ffects 
of the energetically-determined front velocity (4.8). In 
order for the two expressions to be in complete agree- 
ment, the nonlocal contribution of the contour integral 



in (4.8) must then be equivalent to the contribution from 
the inhibition term v^{x = Q) in (4.2S). 



The nonlocal character of the inhibitor field on the pat- 
tern is expressed mathematically by the Green 's function 
solution of the outer inhibitor equation (4.14) 



z;(x) = --^ / dx'Xo(lx-x'l) , 



(4.29) 



where the integration is done only over black areas. Us- 
ing the PDE for Kq and Green's theorem, the above so- 
lution can be rewritten in terms of a line integral over 
the interfacial boundary r. This results in the arclength 
integral 



v{r) 



r r 



-Xi(lr-r'l), (4.30) 



where Ki is introduced via the Bessel identity Ki = 
—Kq. It should be noted that a Plemelj-type argument 
Q is required to evaluate this integral since v° must 
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be evaluated for a point on the boundary r. Here how- 
ever, the usual principal value integral is unnecessary be- 
cau se th e singularity in the integrand is removable (see 
Eq. U). 



Substitution of the contour integral (4.30) into the 
asymptotic front speed ( 4.28 ) confirms the energetic re- 
sult ( |4.4| ). In terms of v^{x — Q), the global energy of a 
black pattern may be written as 



£ = -fL 



w°(x) 



2AF 



(4.31) 



a result that is not immediately obvious from the asymp- 
totic derivation. Written in this way, with the area A 
written as an integral, we see immediately from the fact 
that u° — 1 < that the integrand loses positive definite- 
ness for p > 2AF, which is precisely the criterion for the 
stability of stripes (cf. Eq. 3.11 ) 



As formulated here, the front velocity formulas (4.2S 



and (4.4) both involve the evaluation of the outer in- 
hibitor field, either by boundary inte gration (4.3), or 
direct solution of the bulk PDE (4.14). The boundary 



approach offers a compact, intrinsic description of com- 
plex contours; however, the bulk PDE approach can be 
calculationally advantageous in simple geometries where 
natural eigenfunctions can be constructed. An example 
of the latter is the analysis of a disc-shaped domain. 

Localized Disk Solution: For a black circular spot of 
radius R in an infinite white domain, the outer inhibitor 



v'^ has a radially-symmetric solution of (4.14) in terms of 
modified Bessel functions 



0(._ (l-RK^{R)Io{r) r<R 
^ ^ ' ~ \ RIi{R)Ko{r) R<r ' 



(4.32) 



where continuity of and is satisfied atr = R. Direct 
substitution of v'^{R) into (4.28) gives a nonlinear ODE 



for the purely radial evolution of the contour 

Rt - -6V2dI -\ —- + AF 
[6 V 2 i? 



+P 



Rh{R)Ka{R)-- 



(4.33) 



The dynamics (4.33) shows that the condition for the 
existence and radial stability of equilibrium disk solu- 
tions depends only on two 0(1) parameters: the poten- 
tial difference and the inhibitor coupling, both scaled on 
the activator diffusion length. Equilibrium radii R* then 
satisfy the transcendental condition 



(a) two equilibria, only larger radius stable; 

(b) one equilibrium, unstable; 

(c) no equilibria, radial contraction only; 

(d) three equilibria, intermediate disk stable. 

Apart from the anomalous region (d), note that black 
localized disk patterns can only be stabilized when AF > 
and white is energetically-preferred. This behavior is 
the analog of the one-dimensional stripe stabilization and 
illustrates the process of domain localization. Although 
stable disks can be sustained in the small region (d) for 
AF < 0, since this region is so small in parameter space 
and the attracting basin for the stable radius relatively 
narrow, we believe these states to be unimportant, at 
least, in the context of this particular formulation of the 
FitzHugh-Nagumo system. 

It is also important to note that in the absence of in- 
hibitor coupling, no stable disks are possible. The case 
when AF — (and p = 0) corresponds to the well-known 
shrinkage by curvature that leads to a singular collapse 
in finite time 



R{t) - ^/2D{tn - t) 



(4.35) 



Note that the PDE evolutions shown in Figs. ^ and 
^ both have parameter values that lie within region 
(a). The dynamical evolution in ( 4.34 ) is restricted 
to radially-symmetric geometries; the possibility for 
azimuthal instabilities, which ultimately generate the 
labyrinthine patterns, is addressed in the following sec- 
tion. 

Just as we obtained the energy of a domain in one di- 
mension and deduced the existence of an energy barrier 
to its creation, we may compute the ener gy of a circular 
domain of radius R on the basis of Eq. 4.4 . The computa- 
tion of the nonlocal contribution may be done two ways, 
via the direct calculation of the self-induction integral, or 
by finding the inhibitor field v(r) associated with a cir- 
cular activator pulse. Collecting all of the contributions, 
we obtain the energy of a circle of radius R, 

A£{R) = 2ttR-/ + ttR^AF - irpR^ Ki{R)Ii{R) , (4.36) 



Using the simplest estimate of 7 from Eq. 3.35| , Figure 
shows for se veral different values of p the function A£{R) 
in Eq. i.Zt , illustrating the presence of a local minimum 
at a finite value of R for sufficiently large p. 



V. INSTABILITIES OF CHEMICAL FRONTS 



A. Mechanism of the transverse front instability 



6V2i?* 



AF + p 



R*h{R*)Ko{R*) - 



1 



= (4.34) 



which, upon numerical root tracking, identifies four dis- 
tinct equilibrium scenarios. These are illustrated in Fig- 
ure O: 



The last of the three important aspects of the reaction- 
diffusion dynamics outlined in the introduction concerns 
the instabilities of chemical fronts. At the level of linear 
stability analysis, we will see that the mechanism of this 
instability is rather similar to that of the MuUins-Sekerka 
instability of a liquid-solid interface 1 41 1 . In Figure |ll| we 
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show schematically the level curves of activator and in- 
hibitor concentration near a modulated interface of the 
inhibitor. The excess of inhibitor in the concave regions 
of the front has the tendency to push those concavities 
further back. Conversely, the portions of the front that 
are convex outwards are locally depleted in inhibitor rel- 
ative to the flat interface, and will expand further out- 
ward. Both of these effects increase the length of the 
interface and its curvature, and will thus be resisted by 
line tension. The two effects will balance at a character- 
istic length scale to produce an instability. 

The instability may be understood in its simplest form 
at the level of the energetics of circular domains; con- 
sider again the expression in Eq.( [4.36| ) and examine the 
local limit 1/R <C 1 In this limit, we use the asymp- 

totics Ki{z)Ii{z) ~ (l/2z)(l - 3/8z2 h ) to simplify 

the energy 



AS = 2ttR [l-^P 



ttR^U 



'inp 



(5.1) 



We see two consequences of the inhibitor coupling. First 
there is an effective line tension 7 — p/4 which may be 
driven negative for sufficiently large activator-inhibitor 
coupling. Clearly a negative line tension will lead to pro- 
liferation of the interface. A negative line tension is a 
common feature of fingering instabilities arising from the 
competition between Young-Laplace and Biot-Savart in- 
teractions | |20| , |39| , p^ , p3[ . Second, at higher-order in 1/R 
there appears a stabilizing contribution which below we 
show is like that for the bending energy of an elastic rod. 
If the effective tension is negative, then this last con- 
tribution will stabilize the interface against small-scale 
disturbances and produce a finite-wavelength instability. 



B. Approximate Local Contour Dynamics 

The heuristic notion of an effective line tension may 
be made more systematic by recasting the nonlocal con- 
tour dynamics in an approximate local (but still intrin- 
sic) form through an expansion in powers of the pre- 
sumed small curvature /t, and its arclength derivatives. 
The method of this expansion is similar in spirit to that 
used recently [Q to study screened electrostatic poten- 
tials near surfaces of arbitrary geometry. Indeed, the 
inhibitor field obeys the same modified Helmholtz oper- 
ator as appears in the Debye-Hiickel theory of colloidal 
interactions [Q, which leads to the screened Coulomb 
interaction between elementary charges. 

The curvature expansion may be constructed directly 
from the equation of motion, or first at the level of the en- 
ergy functional and then carried through the variational 
principle. Adopting the latter method, we first expand 
the scala r pr oduct of tangent vectors in the self-induction 
integral (L4) as a power series in A = s' — s. 



t(s) -iis') ~ 1 



(5,2) 



and the distance function 



|r(s')-r(s)|c.|A|--A2|A|«:^ 



(5.3) 



Substituting into the self-induction integrals and extend- 
ing the limits of integration over s' to ±00 (with expo- 
nential accuracy), we obtain 



(5.4) 



ds fds t ■ t Ko ~ n fds [1 n (s) 



and the approximate local energy functional 

A£[r] HA + + (j)dsK^ , (5.5) 

where the effective line tension is as in Eq. |5.1|: 7 = 
7 — {1/4) p. The term proportional to in Eq. p3.5| is the 
bending energy of an elastic line , and its coefficient is 
positive if the activator-inhibitor coupling constant p is 
positive. Under that condition, it prevents the interface 
from bending on arbitrarily fine scales. 

Using the approximate local energy functional ( ^.5| ), 
the normal velocity of the interface is 



7 . 

— n • r* 

D 



-AF 



+ — p 1 



1 . 

— K' 

2 



(5.6) 



Apart from the pressure term, this is the planar version 
of the "curve-straightening equation" |0 that is equiv- 
alent to the Rouse dynamics of the worm-like model of 
elastic polymers |Q . It has the form of a Landau expan- 
sion in powers of the curvature and its derivatives, and 
in that sense is similar to so-called "geometric" models 
of crystal growth [Q. Unlike those mode ls, h owever, the 
coefficients of the various terms in Eq. 5.6 are not all 
independent. Thus, for instance, the terms Kss and 
must have coefficients with ratio 2 by virtue of the vari- 
ational principle applied to the energy functional ( ^.5[) . 



C. Linear stability for fronts and stripes 

The most fundamental instability of chemical fronts in 
the present model is that of an infinite straight interface 
bounding the states u — and u = 1. To compute the 
spectrum of growth rates of perturbations to this front, 
we parametrize the interface as 



r(s) = xe^ + C{x, t)ey , 



(5.7) 



If we let C{x,t) — Ckit) cos{kx), the linearization of the 
Biot-Savart integral is 



rfs'R X tKi ~ Q cos{kx) X 
dy 



y 



(fcysin(2fcy)-2)Xi(y) , (5.8) 
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Use of the identity {\/y)Ki{y) = -Kpiy) - K[{y) and 
an integration by parts transforms (5^) into a standard 
integral. Combining these results with the linearized cur- 
vature, the growth rate a{k) — {dtCk)ICk for the single 
front instability is 




1 



(5.9) 



(5.10) 



For wavelengths that are small relative to the 0(1) in- 
hibitor screening length, we obtain 





1 


{- 





(5.11) 



showing, as disc ussed in the previous section (see in par- 
ticular Eq. 5.6), that the instability arises from a neg- 
ative effective line tension, stabilized by an elastic-type 
term in k^. This result also shows that at the onset of 
the instability the critical wavelength is infinite. As in 



the analysis of steady states for disk solutions, Eq. [1.34 
we shall find it convenient to work in a rescaled param- 
eter space with coordinates p and f as in ( |3l3| ). Setting 
a{k) = and da{k)/dk^ = we obtain the critical cou- 
pling constant for the instability. 



Pf 



47 _ 72 



(5.12) 



Note that this is independent of the parameter r, and 
hence also of the energy difference AF, and so is equally 
valid as a stability criterion for a uniformly moving 
straight front. The critical value of p could also be ob- 
tained directly from the vanishing of the effective line ten- 
sion m ( PI ). For p > Pf, the most unstable wavevector 
is 



k* = 



2/3 



-,1/2 



(5.13) 



Figure |l2| shows in the rescaled parameter space the lo- 
cation of the instability of a single front, as well as the 
boundary of stable stripes, defined by 



1 



P - 



(5.14) 



Beyond the behavior of a single chemical front, it is 
natural to consider the stability of stripes. Perturbations 
to the shape of a stripe may be decomposed into those 
with odd and even parity under reflection through the 
stripe midline, termed sinuous (S) and varicose (V), re- 
spectively. In these two situations, the interfaces bound- 
ing a stripe of width 2Q* are parametrized as follows: 



r±{x) =xe^ + {±Q* +({x,t))ey (S) 



(5.15a) 



r±{x)=xe, + {±Q*±C{x,t))ey (V). (5.15b) 

A straightforward calculation yields the growth rates of 
these two modes: 



(js(fc) = crf(fc) - 3V2l3pe-^^' 

1 



X 1 



/(fc) 



-2Q*(/(fc)-l) 




-2Q*(/(fc)-l) 



(5.16a) 



(5.16b) 



In each case, the growth rates are less than that of the 
single front. The sinuous mode grows faster than the 
varicose, whose damping at fc = reflects the local sta- 
bility of the preferred stripe width. This is consistent 
with what is known in other contexts, for instance in the 
buckhng instability of magnetic stripes Figure |l^ 

shows these growth rates for typical values of the param- 
eters. The greater stability of a stripe relative to a single 
front may also be revealed by the onset of its sinuous in- 
stability at higher values of p. Expanding the growth rate 
( ^.16a ) to order /c^ as in Eq. (5.11), we obtain marginal 
stability conditions by setting the coefficient of fc^ to zero, 
yielding the transcendental equation 



V2 r f 
-^'V-Tp 



l + ln(^ 



(5.17) 



A numerical solution to this is shown in Fig. |l^, where 
we see that this stability boundary properly merges with 
that of a single front as r — > 1/2"*", and the equilibrium 
stripe width Q* diverges. Note that Eq. 5.17 is precisely 
the condition that the structural energy per unit length 
of the stripe (|3.12|) vanish. That these two conditions 



are related is seen by the fact that the coefficient of k in 
the stability analysis is an effective line tension or energy 
per unit length. Beyond this stability boundary, the pro- 
liferation of stripes (straight or buckled) is energetically 
preferred. 



D. Azimuthal Instability of Disks 



The evolution equation (4.32) that establishes the size 



of equilibrium disk solutions (4.34) also determines their 
nonlinear stability to radial perturbations. For azimuthal 
disturbances, the growth exponents are obtained by con- 
sidering a disk whose radius develops a small sinusoidal 
variation 

r = R* +(^{t)smn0 . (5.18) 
The linearized curvature of the interface as calculated 



using (4.24) is 
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R* ^ ' R*^ 



(5.19) 



where it is noted that the n = ±1 perturbations are 
equivalent to simple translations of the disc. In this near- 
circular geometry, it is advantageous to use the bulk ap- 
proach to solve for the inhibitor field since the leading- 
order correction to the unperturbed inhibitor field in- 
volves only a single eigenmode. To 0{C}, the inhibitor 
field v^{r) is 



CR*K*I„{r)smn0 r<R 
(R*I*Kr,{r)simie R<r 



(5.20) 



where u°(r) is given in (|3|), K* = Kn{R*) and /* = 
In{R*), and where the continuity conditions must now 
be imposed on the perturbed interface. This gives the 
corrected value of the inhibitor on the interface 

u°(r) = v°iR*) - (R* [Kin - lnK\ sin"0 , (5.21) 

where the correction also vanishes for the n — ±1 modes. 

Collectin g th e 0(C) corrections from the front speed 
dynamics ( 4.28| ) yields the linearized growth exponent 
for the n*'^-mode perturbation 



7- 



R* 



pR* [Kill - i:ki] 



(5.22) 



The neutral curves for the lower mode instabilities acting 
on the disc solutions of region (a) are shown in Figure |^. 
Coincidentally, note that the conver genc e point for these 
curves is identical to the critical pc ( |5.12 ) which is where 
R* is naturally large. The fact that the instabilities of 
disks and straight fronts derives from a common mecha- 
nism is demonstrated in the limit of large disks yet finite 
wavenumber k = n/R*, yielding 



lim 

n,R* — 



f^(fe=n/-R*) = ^/(^) 



(5.23) 



VI. NUMERICAL STUDIES 
A. Simulation of PDEs 

In this section we presen t nu merical studies of the 
reaction-diffusion dynamics (2J) to verify the stability 
results discussed in the previous section, within the con- 
text of the parameter space of Fig. [ij. In the subsequent 
section we investigate the contour dynamics for compar- 
ison. These numerical studies were performed using a 
pseudospectral algorithm in a two-dimensional periodic 
domain. The algorithm is outlined in Appendix A. In 
showing the time evolution it is convenient to adopt a 
rescaled time r = \j2D{r — l/2)t 

The first phenomenon illustrated is the two-dimen- 
sional version of the approach in one dimension of two 



chemical fronts illustrated in Figure |5[ This is shown 
in Fig. 14 starting from an initial condition in which 
the two fronts are far apart and given a modulation by 
a random collection of Fourier modes. To indicate the 



symmetry of the underlying PDEs given in Eq. 2.5, the 
parameters chosen {D — 0.01, p = 0.2, r — 0.4) are such 
that the state m = 1 is the more stable and invades the 
regions with u — 0. The straight front is linearly sta- 
ble, so in addition to the net motion of the fronts toward 
each other, these initial modulations relax. As in the one- 
dimensional studies, for these parameter values the fronts 
do not cross, coming to rest at a finite distance set by the 
detuning r — 1/2 and the inhibitor coupling p. Similar 
behavior has been seen in the iodide-ferrocyanide-sulfate 
reaction [0, as well as in simulations of the Grey-Scott 
model ||]. 

The growth of a labyrinthine structure from a com- 
pact initial condition was illustrated in Fig. |^. The pa- 
rameter values for this computation, f = 0.2, p — 1.5, 
are in the region of rescaled parameter space in which 
the simple disk solution is unstable to many azimuthal 
modes. The black regions in the figure are those within 
which u > 1/2. We see in the early stages of the evolu- 
tion the growth of fingers of a well-defined width ~ this 
is the behavior discussed by Ohta et al. [||. Several of 
the fingers undergo tip-splitting. Their mutual repulsion 
leads eventually to a spacc-filling labyrinth which appar- 
ently converges to a steady state. This convergence is 
clearly a consequence of the front self-avoidance in the pe- 
riodic computational domain. The interactions between 
the fronts have been such as to create a rather uniform 
width to the fingers of u = 1, as well as to the intervening 
regions with u = 0. Besides their similarity to the exper- 
imental patterns of Lee, et al. ||l|,|^, the phenomenology 
of this pattern formation has a very strong resemblance 
to that seen in magnetic fiuidsjl^, superconductors pi] , 
as well as thin garnet films |2j ]. 

By changing the coefficients to f = 1, p = 1, we en- 
ter the regime in which the localized disk is stable, and 
a branched structure may relax to it without fissioning. 
This is shown in Fig. ^, in which the starting field con- 
figuration is panel (c) in Fig. ||. Again, this shape relax- 
ation is like that seen in the iodide-ferrocyanide-sulfite 
reaction, as well as that observed in magnetic fiuids in 
Hele-Shaw fiow when the applied magnetic field is re- 
moved and surface tension returns a fingered structure 
to the circular ground state pct |. 

The front interactions responsible for labyrinthine fin- 
gering instabilities of a single domain naturally appear 
in the interactions of multiple domains, and can lead to 
patterns composed of disconnected but highly interdigi- 
tated regions. Figure ^ shows contour plots of the evo- 
lution starting from an initial condition with two small, 
nearly circular domains. Despite the complex fingering 
instabilities, the domains retain their integrity and do 
not merge. This multiple-domain problem is rather sim- 
ilar to that observed in the inter mediat e state of Type-I 
superconductors (see also Section VIII ) . 
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An important aspect of the evolution of compact do- 
mains concerns the possibihty of domain fission and fu- 
sion. As is quite common with interfacial treatments of 
pattern formation, the contour dynamics is not asymp- 
totically valid when regions of the interface approach each 
other near the reconnection point point. The full PDEs 
for the activator-inhibitor pair are however well-defined 
during these events, and both fission and fusion are possi- 
ble if the pressure driving the interfaces together is large 
enough, or conversely the inhibitor-induced repulsion is 
small enough. This has also been observed in experi- 
ments 101. 



B. Simulation of contour dynamics 

Next we turn to numerical studies of the nonlocal con- 
tour dynamics. The numerical method for this, outlined 
in Appendix A, has been employed elsewhere for 
analogous problems in nonlocal interface motion. It is 
based on a pseudospectral treatment of the tangent-angle 
representation of the dynamics. 

The most elementary instability of a localized state is 
the elliptical one, illustrated in Fig. Here, the initial 



condition is a solution of the equilibrium condition (4.34) 
for a localized disc, perturbed with a small amplitude of 
the n — 2 mode. After a short transient, the domain 
enters an era of linear elongation as illustrated in the 
lower panel of the figure. This linear growth is a r eflection 
of the negative structural energy density ( |3.12D for the 
parameters f = 0.212, p — 0.796. 

Increasing the value of the inhibitor coupling to p = 
0.90 renders the localized state unstable to the mode 
n = 3 (as well as to the elliptical mode). The bulbous tips 
that form during the elliptical instability have a nonzero 
projection onto the n = 3 mode, leading to tip-splitting 
events. This is shown in Fig. |l^. In this regime of 
parameter space each time a new tip is formed it is sus- 
ceptible to tip-splitting, leading to a cascading process 
and a proliferation of three-fold coordinates vertices. 

For p larger still, we may investigate one of the com- 
mon features of labyrinthine structures seen in dipolar 
systems as well as in chemical reactions: the appearance 
of three-fold nodes to the exclusion of all others of higher 
coordination. For certain kinds of optimization problems 
involving the minimization of interface length under the 
constraint of fixed endpoints it is known that three-fold 
nodes are the only stable vertices , but no such result 
is known in the present context. The numerical simu- 
lations strongly suggest, however, that higher-order ver- 
tices are dynamically unstable. Figure |l^ shows the evo- 
lution of an initial condition consisting of a localized disk 
modulated by four-fold perturbations of several per cent 
and a two-fold distortion one fiftieth as large. A cross- 
like vertex forms quickly, but is unstable to the elliptical 
perturbation, splitting into two three-fold vertices. 
Having seen in isolation the elementary processes un- 



derlying labyrinth formation, stripe formation and pro- 
liferation, tip-splitting, and vertex reduction, we show in 
Fig. |l^ the appearance of a labyrinthine pattern from a 
compact initial condition. The characteristic feature of 
interface repulsion and the appearance of well-defined fin- 
ger widths is readily apparent. Unlike the simulations of 
the PDEs, the contour representation does not naturally 
build in periodic boundary conditions, so the long-time 
evolution of the two will differ significantly. Changing 
the coefficients so that the localized disk is both radially 
and azimuthally stable, we see in Fig. ^ the relaxation 
to a compact state of a branched initial condition taken 
from near the the end of the evolution in Fig. p9l 



VII. BEYOND THE SLAVING LIMIT 

As we have discussed so far, the limit e — > renders the 
dynamics an overdamped gradient flow, associated with 
a Lagrangian variational principle in which the kinetic 
energy is neglected and the dissipation function is local. 
As the recent work of Hagberg and Meron [p6|j2^ ] has 
emphasized, the nature of the chemical fronts between 
the two metastable values of the activator may be quite 
different in the two extremes of slow- and fast-inhibition. 
From Eq. 4.25, the timescale of the front dynamics is 



Ut — 0{D), so that in actuality, the fast-inhibitor limit 
only requires that e 13 <C 1 - so that the value of e can be 
quite large. Extending this contour dynamic approach 
beyond this fast-inhibitor limit represents a significant 
challenge. Nevertheless, an intuition for how oscillatory 
behavior might arise with finite e is suggested by argu- 
ments using contour energetics. 



For e 7^ we may solve the inhibitor dynamics (2.4b) 
for V in terms of u, 



u(x,i) 



J dt' J dx'G{x -x',t- t')u{x', t') 



where for i > the Green's function is 



G(x,t) 



exp 



ex 
IF 



(7.1) 



(7.2) 



Clearly, e is the natural time scale for the decay of G, and 
if we take the liberty of expanding the integral for slowly- 
varying u we obtain u(x', t') ~ u(x', t) + {t' — t)ut{x' , t) + 
{l/2){t' ~tfuu{yi' ,t) + - ■ •, substitute into (O) and per- 
form the time integrations. Each power of t' — t in the 
expansion will contribute a power of e. Up to order 
we obtain 



w(x,t) 



J dx'|go(x ~x')u{x',t) + egi(x - x')Mt(x', 



+e2g2(x-x')u„(x',t) 



(7.3) 



where Go{:x.) = {1/2t:)Kq (x) and the remaining functions 
G, are 
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^i(x) = -^xKi (x) , g2(x) = -r^x^K2 (x) . (7.4) 

47r lOTT 

Note the very important feature that the order e kernel 
Qi is negative, while the second-order kernel is positive. 
Up to quadratic order in e the nonlocal activator dynam- 
ics may be written in the following form 

/SS 
(ix'^2(x - x.')utt{x.' ,t) + — = 
du 

-ut-ep y"dx'gi(x-x')ut(x',t) . (7.5) 



This conforms to the variational principle in Eq. (3.7) 
having the form 



di 



5T 

Sut 



5£_ 

6u 



sn 

Sut 



(7.6) 



but now with a finite kinetic energy functional 

r^ie^p Jdx y"dx'ut(x,i)52(x~x')wt(x',t) , (7.7) 
and a nonlocal contribution to the dissipation function 

+^ JdxJ dx'Mt(x,t)^i(x-x')ut(x',i) . (7.8) 

The positivity of Q2 renders the "mass" of the field dy- 
namics positive, while there is a competition between the 
signs of the local and nonlocal parts of the dissipation 
function. 

A heuristic derivation of the leading-order changes to 
the contour dynamics from the inclusion of nonzero e pro- 
ceedsas in Section IV. Starting from the dyn amics for u 



in (7^) or equivalently the functionals in (7/7) and (7^ 



we recognize that the time derivative Ut and acceleration 
Utt are localized at the domain boundary. Using a cor- 
respondence like that in Eq. 3.9, we obtain the contour 



representation of the kinetic energy, 

T ~ ^ds jds'-a ■ vt g2{v - r') n • v[ 



(7.9) 



and dissipation function 



7^ - ^ i>ds{h-vtf 



<j)ds tj)ds'h ■ VtQi{Y — r')n • 



(7.10) 



This kind of nonlocal dissipation function occurs as well 
in the analysis of certain models of solidification |^,^ , 
and is conceptually similar to role of the Oseen tensor in 
polymer dynamics [ p3| . 

The interpretation of this dynamics as a damped me- 
chanical system is complicated by the nonlocal nature of 



the kinetic energy and dissipation functionals. With the 
function Qi < 0, the overall rate of dissipation is positive 
for small e but may become negative for e sufficiently 
large, say of order 1/p. One expects qualitatively new 
behavior in this limit. Moreover, the inclusion of inertial 
terms in the dynamics makes the analysis of traveling- 
wave states qualitatively different from that in the over- 
damped limit. Specifically, an ansatz of the form u{x—ct) 
yields a cubic equation for the speed c, with three pos- 
sible roots, in contrast to the unique speed found in the 
limit e — > 0. These two roots most likely represent those 
associated with the nonequilibrium Ising-Bloch front bi- 
furcation discussed by Hagberg and Meron and are 
related to the oscillatory instabilities discussed by Ohta, 
et al. ^ for the case of simple geometries. 

Developing a theory for the asymptotic stability of the 
labyrinthine patterns is one current effort [|4| . It is hoped 
that results in this direction will lead to a quantitative 
connection between the problem of finite inhibitor diffu- 
sion (e 7^ 0) and the onset of time-oscillatory behavior. 



VIII. NONLOCAL CONTOUR DYNAMICS IN 
OTHER SYSTEMS 

We have focused here on the fast-inhibitor limit of a 
reaction-diffusion system and demonstrated that its be- 
havior is well-described by a nonlocal contour dynamics 
model. In this limit, where the inhibitor is slaved to the 
activator, the dynamics is a gradient flow with an energy 
functional of the form 

£[y] =nA + jL-^p jds jds'i ■ t'$ {R/h) , (8.1) 

with R = |r(s) — r(s')|, and a normal velocity U propor- 
tional to the force obtained variationally as — n • S£/Sr, 
with 



[/ = -n-7K 



J f ds'R(s, s') X t'$' {R/h) . (8.2) 



In this section we make specific the connections between 
this form of dynamics and those found in several other 
quite distinct physical systems. These are (i) Type-I su- 
perconductors in the intermediate state IpqL (ii) mag- 
netic fluids in Hele-Shaw flow (l|j2^,||,|||6l and (iii) 
Langmuir monolayers of dipolar molecules 1 43|^ 5^ . While 
the functions <I> differ from one system to the next, the 
common feature we find in each is a positive bare line 
tension 7, and a repulsive interaction between antiparal- 
lel tangent vectors, with p > 0. Interfacial instabilities 
leading to fingered structures are then seen to derive from 
a negative effective line tension arising from the nonlo- 
cal contribution. It is this same nonlocal coupling of the 
interface that leads to the self-avoiding nature of the pat- 
tern formation beyond the linear instability. We discuss 
each of the three systems below, pointing out the dif- 
ferent physical origins of the nonlocal coupling, and the 
different constraints on the interface motion. 
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Type-I Superconductors: The intermediate state of 
Type-I superconductors occurs when a thin slab of the 
material below its zero-field transition temperature is 
placed in a magnetic field normal to its surface [^ . 
Rather than exhibit a complete Meissner effect, the de- 
magnetizing effects arising from the sample geometry 
lead instead to the penetration of the flux through the 
sample in an intricate arrangement of flux domains, each 
of which is fingered and often branched. The shapes of 
these flux domains arise from a competition between the 
positive superconductor-normal surface energy and the 
interactions between the Meissner currents which circu- 
late at the flux domain boundaries. As argued many 
years ago [|8 59|, this self-induction interaction retains 
its long-range form as in free space, since the electro- 
magnetic fields in the vacuum above and below the slab 
are unscreened. In the simplest treatment of these in- 
teractions ll^, they are considered identical to those of 
current loops in completely free space. Upon averaging 
the standard Coulombic self-interaction between elemen- 
tary current segments over the thickness h of the slab, 
one finds a tangent-vector coupling for a single interface 



$(z) 



In 



z-' + il + z-^V^^' 



-<f'iz)=[l + z-'f'-l, 



+ z-[l + z'f' , 
(8.3) 



where the characteristic length scale is h. While retaining 
the Coulombic form ^{z) ~ 1/z for z ^ oo, <I> has only a 
logarithmic singularity at the origin, much like the Bessel 
function Kq in the reaction-diffusion problem considered 
here. 

Interface motion in superconducting systems does not 
conform to the simple local-dissipation model discussed 
in Section IV, but rather reflects the diffusion of the mag - 
netic flux in the normal phase Thus, the force ( |8.2| ) 
becomes a boundary condition for the diffusion equation 
obeyed by the field, rather than determining the velocity 
directly |38| . In simple models of this pattern forma- 
tion p5| , in which the local dissipation model is invoked, 
this conservation law leads to the generalization of the 
product HA to the derivative of a bulk free energy den- 
sity Sbuik (A) which arises from the competition between 
the field energy of the external magnetic field and the 
condensation free energy of the superconducting state. 
The conservation of magnetic fiux is an important global 
constraint in this problem, and leads to an equilibrium 
area fraction of the domains that is determined primar- 
ily by the minimization of this energy; the shapes of the 
individual interfaces arise from the competing line ten- 
sion and Biot-Savart interactions. Unlike the reaction- 
diffusion problem, the intermediate state properties are 
determined fundamentally by many-domain interactions. 

Magnetic Fluids: A second physical system conform- 
ing to the energetics in Eqs . |8.1| and S.2 is that of thin 
domains of magnetic fluids ]17|| in the geometry of Hele- 
Shaw flow. There, the domain is trapped between two 
glass plates spaced a distance h apart, with the remainder 



of the gap filled by water. A magnetic field applied nor- 
mal to the plates aligns the microscopic magnetic domain 
in suspension, creating an approximately uniformly mag- 
netized droplet. By the correspondence between mag- 
netization and current-loops, the field energy associated 
with the domain is again represented in terms of a self- 
induction interaction, with the functions $ and ^! as in 
(^.3| ) by virtue of the slab geometry. The amplitude p is 
now proportional to M^/i, where M is the magnetization 
of the domain, and the line tension 7 = ha, with a the 
ferrofluid-water surface tension. 

In this hydrodynamic problem, the dynamics is well- 
approximated by Darcy's law v —(h^/12ri)'VP, where 
V is the z-averaged in-plane fluid velocity, 77 the fluid vis- 
cosity, and P a generalized pressure including magnetic 
contributions |^^,^ . With the constraint of fluid incom- 
pres sibility, the pressure fleld is harmonic, with the force 
( ^.2| ) again acting as a boundary condition on P. The 
fluid incompressibility leads directly to the conservation 
of the area enclosed by the boundary. 

Langmuir Monolayers: The flnal class of systems gov- 
erned by these energetics includes amphiphilic (Lang- 
muir) monolayers at the air-water interface. These con- 
sist of single- or multi-component monomolecular films 
of surfactants in which the mean lateral density is regu- 
lated externally, allowing for a study of phases and phase 
transformations. Under suitable conditions, domains of a 
high-density phase appear in a background of lower den- 
sity, and may be visualized by the differential fluorescence 
of a dye incorporated into the layer. One observes various 
shape instabilities of these domains as conditions such as 
temperature and pressure are varied fl^ . These are be- 
lieved to arise from the competing effects of line tension 
at the domain boundary and long-range electric dipole 
interactions between the molecules which are oriented by 
the constraint of packing. The energetics of these inter- 
actions may be described by taking the "ultra-thin" 
limit of the magnetic formulation (B.3) above, with the 
cutoff h being a molecular length, and the amplitude p 
proportional to the square of the dipole moment density 
p. In the limit of small h, the self-induction and Biot- 
Savart terms have the form associated with infinitesimal 
current-carrying wires |57[ , 



a>(z) 



-$'(z) 



1 



(8.4) 



A cutoff procedure must be implemented on these func- 
tions to treat the divergences which occur when s — 
s' — > 0. As in the superconductor problem, in the sim- 
plest model [|43| 11 is a Lagrange multiplier conjugate 
to the area. Experiments have shown a variety of 
fingered domain shapes consistent with the boundary 
model. More recent work has focused on the hydrody- 
namics of monolayer domains coupled to the viscous sub- 
fiuid (6||6|. 

The nonlocal energy functional £[u] in Eq. ( 2.18| ) 
is known also to be relevant for physical systems quite 
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distinct from those with dipolar interactions, appear- 
ing for instance in models of microphase separation in 
block copolymers [Q. There the nonlocal coupling is 
long-ranged, reflecting the connectivity of the polymers. 
Labyrinthine patterns occur there as well (see also the re- 
view in Ref. [P3| ), but are not necessarily confined to two 
dimensions. It has also been remarked [6^ ] that this non- 
local interface coupling may be related to the diffusion 
of impurities and/or latent heat in solidification, which 
leads to some degree of interface self- avoidance seen in 
the development of dendrites. Such behavior can not 
be captured by purely local geometric or boundary- 
layer |6^] models. Finally, observe that the Biot-Savart 
coup ling a lso appears in the contour dynamics formula- 
tion [ |67t|6^j69tl of vortex patch motion in two-dimensional 
ideal fluids. Rather than exhibiting strongly overdamped 
dynamics, these are of course Hamiltonian systems. 



IX. CONCLUSIONS 

Finally, we discuss briefly three important open issues 
regarding the dynamics of labyrinthine pattern forma- 
tion: generalizations to higher spatial dimension, further 
elucidation of variational principles, and derivation from 
microscopic chemical kinetics. 

We have seen that the phenomenon of lateral inhibi- 
tion necessarily involves interactions between segments 
of chemical fronts that are potentially far apart in ar- 
clength, yet close in space. It i s pr ecisely this nonlocal- 
ity which enters the energetics (3.1) and dynamics (B.2). 
Whereas in the magnetic systems the angular part of 
the current-current interactions is naturally written as 
t(s) • t(s') in accord with the existence of currents circu- 
lating in the direction t, there is no such circulation in 
the reaction-diffusion problem. But of course, as seen in 
the derivation (4.3), this scalar product may equally well 
be written in terms of the normal vectors as n(s) • n(s'). 
This formulation makes it clear that the central issue 
is whether fronts oppose one another with a region of 
low activator concentration in between. The normal vec- 
tor representation allows a straightforward generalization 
to the interaction of two-dimensional surfaces, for which 
there is no natural or unique assignment of tangent vec- 
tors. We suggest that certain three-dimensional patterns 
may be profitably studied by models embodying this non- 
local interaction. In addition to the block copolymer sys- 
tems mentioned earlier, other candidates include highly 
convoluted structures such brain coral, which displays 
labyrinthine structures with features such as three-fold 
coordinated nodes like those seen in the present work. 
One may imagine that these arise from the interplay be- 
tween growth of individual members of the colony and 
the competition for nutrients. 



As discussed in Section VII, a heuristically-derived 



variational principle much like that of a damped mechan- 
ical system. An issue of some significance is the extent to 
which the contour dynamics approach may be more rigor- 
ously extended to incorporate the front bifurcation that 
ultimately occurs for e sufficiently large. Related issues 
concern the connection between such a description and 
spiral-wave behavior, as well as the nature of the front 
bifurcation for surfaces moving in three dimensions. 

In light of the present derivation of the contour dy- 
namics from the FitzHugh-Nagumo model, it remains of 
great interest to investigate as well whether the particu- 
lar chemical kinetics pQ] relevant to the experiments of 
Lee, et al. may be recast as an interface dynamics. Re- 
cent work JtH has shown that the those very complex 
kinetics have dynamics on many time scales and may 
be reduced to an effective two-variable model through a 
sequence of slaving approximations. While the form of 
that reduced description is somewhat different than the 
FitzHugh-Nagumo model, the possibility that it shows 
similar, near gradient-flow behavior is an intriguing area 
of investigation. 
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APPENDIX A: NUMERICAL METHODS 

Here are summarized established pseudospectral meth- 
ods ||7|,^ that we have adapted to study both the 
reaction-diffusion dynamics and the contour evolution. 
For the case of a scalar partial differential equation in 
1 + 1 dimensions, 



at 



(Al) 



where £ is a linear operator and A/" is nonlinear. We 
assume that the highest-order spatial derivative appears 
in the linear operator. In Fourier space ( |Al| ) is 



contour dynamics for small deviations from the fast- 
inhibitor limit appears to conform to a rather general 



du{k, t) 
dt 



u}{k)u{k,t) ^Af{k,t). 



(A2) 
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with uj{k) = C{ik). The nonhnear terms are obtained 
pseudospectrally by fast Fourier transforming M in real 
space, N{k,t) = T[N{u{x,t))], and u{x,t) is computed 
by an inverse fast Fourier transformation of u{k, t). Now 
define 



v{k,t) 



- p-'^(fc)* 



u{k,t), 



(A3) 



and multiply (|A2[) by the exponential factor, yielding 



dv{k,t) 
dt 



•J{k)tj^ 



(A4) 



In the simplest Euler method, the solution to Eq. A4 is 



v{k,t + At)-v{k,t) 
At 

which yields 



W (v{k,t)e 



u{k, t + At)= e'^^'^)'^* u{k, t) + AW{k, t) 



(A5) 



(A6) 



The exponentiation of the growth rate uj{k) in ( [Aq ) plays 
a useful role in guaranteeing stability for a diffusive linear 
operator (u;(fc) = —jk"^). The usual stability considera- 
tions would require a time step At such that for 
large values of momentum (near the Brillouin zone edge 
fcmax = 7r/a, with a the lattice spacing in real space) 
the quantity k^At be less than unity. This requires an 
extremely small time step, rendering the calculation pro- 
hibitively slow. Here, even if k^^^At > 1, the calculation 
is stable due to the incorporation of the exact dynamics 
of the linear operator, namely the exponential damping 
at high momentum. 



The generalization of Eq. A6 to a 4th order Runge- 
Kutta method proceeds as follows. Define 



ii/2(fc) ^e'^W'^*/^, L(fc) 



and the intermediate results 



Ui{k) = L,/^{k) 



-AW{i 



U2{k) ^ Li/2ik)u+ -AWi 



(A7) 

(A8) 

(A9) 
(AlO) 



[/3(fc) =i(fc)u + AtLi/2(fc)AA2 , 

with J\fj — Af^Uj^ . Then the time-stepping routine 
analogous to (A6) is 



u(t + At) = L(k)u + -AtL{kW + -AtLi/^{k)Ni 
6 3 

+iAiLi/2(fc)A^2 + ^AtA^3 . (All) 
3 ' 6 

To generalize this method to situations with n coupled 
variables Ui (i = l,2,...,n), we write the equation of 
motion in Fourier space in vectorial form 



Uf = O • u + TV", 



(A12) 



with J7 a matrix of wavevector-dependent growth rates in 
Fourier space and J\f a vector of nonlinear terms obtained 
pseudospectrally. Now suppose that det Jl — tjl = has 
as solutions eigenvalues oji and associated eigenvectors 
ei (i = 1, . . . , n). The matrix T whose columns are the 
components of the eigenvectors defines a linear transfor- 
mation between u and the vector a of expansion coeffi- 
cients in the basis of eigenvectors. That is, u = T • a, 
where T^^ • T = I. The equation of motion (A12) then 
becomes 



(A13) 



The matrix D = T^^ • r2 • T is diagonal. Now let v 
exp(— Dt)a, so 



Vf = e 



(A14) 



It follows that for an Euler method, one need only know 
the matrix 



L = T • e^'^'T 



DAtiTi-l 



for then 



u{k,t + At) = L- u{k,t) + AtA/{k,t) 



(A15) 



(A16) 



whereas in the Run ge-K utta method we obtain the vec- 
torial analog of Eq. 



All 



u(k,t + At) = L • u(fc, t) + -AtL ■ N + i AtLi/2 • Ni 

6 3 ' 

+-A<Li/2-AA2 + -AiA/'3 . 
3 ' 6 

where 

and the intermediate results are 



ui = Li/ 



2 ■ 



u+-AtAf{u) 
'-AtAfi 



U2 = Li/2 • U - 

Z 

ua = L • u + AtLi/2 • 7V'2 



(A17) 
(A18) 

(A19) 

(A20) 
(A21) 



with JVj = N (Uj). 

The fast-inhibitor limit: In the limit e = 0, Eq. ( ^.17 ), 
the single equation of motion for u contains a contribu- 
tion that while nonlocal is nevertheless linear. Since it 
is a convolution, it is local in Fourier space and can be 
incorporated directly into the linear operator C The re- 
sulting transform is 



uj{k) 



-Dk^ 



P 



1 + F ' 



(A22) 



precisely the growth rate of the mode + in Eq. |2.7| . The 
diffusive contribution —Dk^ dominates at large wavevec- 
tor. 



20 



Contour dynamics: In numerical studies of the contour 
dynamics we have employed techniques described else- 
where 1^,^ , summarized briefly here. Starting from an 
equation of motion in the form 



rt = Uh + Wt , 



(A23) 



we study the evolution of the tangent angle |49| 9{s) 
related to the curvature by k(s) = 89/ ds, WOl 



de 

dt 



dU 
'~d7 



kW . 



(A24) 



The spectral method described above requires that we 
utilize a periodic function. A convenient choice is the 
deviation ip of 9 from the linear form 9 = 2tts/L for a 
circle, 



9{a,t) = 27ra + ipia,t) 



(A25) 



A natural choice of gauge is that of "relative arclength," 
for which equally spaced points in the parametrization 
a — s/L remain equally spaced in time. This corresponds 
to a tangential velocity [^,^ 49 



W{a) ^ L{a da' kU - da' nU 
JO Jo 



(A26) 



The highest-order arclength derivative in the ip dynamics 
is then diffusive, 



9-0 
'dt 



7 d'^tp 



+ 



(A27) 



and amenable to the integrating factor method outlined 
above. The dynamical variables of the problem are then 
ip and the contour length L which obeys the simple evo- 
lution 



dsnU . 



(A28) 
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FIG. 1. Chemical pattern formation in the iodide-fer- 
rocyanide-sulfite reaction of Lee, McCormick, Ouyang, and 
Swinney [^]. Low and high pH regions appear respectively as 
white and black by means of a pH indicator. Times proceed 
from upper left to lower right in hours following a perturba- 
tion. Figure courtesy of Lee, et al, and reproduced from Ref 
0. 



FIG. 2. Simulation of the reaction-diffusion PDEs showing 
a compact initial condition undergoing a fingering instabil- 
ity, eventually producing a space-filling labyrinth. Panels are 
contour plots with it < 0.5 shown white, and u > 0.5 shown 
black, at rescaled times of r = 0, 5, 10, 15, 20, 25. The param- 
eters in Eq. are D = 0.01, p = 0.15, and r = 0.52; thus. 



the state u = 1 (black) is less stable than u = (white). 



FIG. 3. Reaction-diffusion simulations as in Fig. g, but 
starting with a branched domain, and with D = 0.01, 
p — 0.10, and r — 0.60. With these parameter values, a 
circular localized state is the stable configuration to which 
the system relaxes. 



FIG. 4. Regions of the r — p parameter space in which Tur- 
ing bifurcations are possible (above solid lines) and forbidden 
(below) within linear stability. Inset shows growth rate of the 
most dangerous mode as a function of wavenumber both in 
the linearly stable and unstable regions, for values indicated 
by the symbols. 



FIG. 5. Space-time portrait of the formation of a localized 
state in one dimension by the shrinkage of a region of it = 1. 
Activator is shown as solid line, inhibitor as dashed line, with 
time increasing downwards. Solid symbols locate the value 
It = 1/2. Parameter values are; D = 0.01, r — 0.55, p = 0.15. 
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FIG. 6. Space-time portrait of the growtli of a localized 
state in one dimension from a small nucleated domain. Acti- 
vator is shown as solid line, inhibitor as dashed line, with time 
increasing downwards. Parameters are: D = 0.01, r = 0.52, 
p = 0.15. 



FIG. 13. Spectrum of growth rates for instabilities of iso- 
lated planar fronts, and the sinuous and varicose instabilities 
of stripes of finite width. Parameters are D = 0.01, p = 0.3, 
and 2(5* = 1.5. The single front is more unstable than the 
finite- width stripe, while the varicose mode is the most stable, 
particularly as fc ^ 0. 



FIG. 7. Front dynamics, (a) Velocity of a single front as a 
function of the activator-inhibitor-coupling, for diffusion con- 
stant D — 0.0001, normalized by the velocity for p = 0. 
Symbols indicate results of numerical simulations of the reac- 
tion-difi^usion equations, while solid line are predictions from 
the asympto tic an alysis with the corrected friction coefficient 
given by Eq. |j.30| . (b) Asymptotics for separation versus time 
for two approaching fronts in one dimension, compared with 
simulation results for D = 0.0001, r = 0.515, and p = 0.02. 
Dashed line is prediction of asymptotics with bare friction 
coefficient, while the dotted line includes order p correction 
to front profile and compares well with the full numerical so- 
lution to the reaction-diffusion equations (solid). Inset (c) 
shows the equilibrium separation of the fronts versus p for 
both simulation (solid symbols) and asymptotics (solid line) 
from Eq. 
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FIG. 8. Illustration of the Biot-Savart interactions of 
nearby chemical fronts, as described by the nonlocal contri- 
bution to the force law in Eq. 4.8 . 



FIG. 9. Regions of existence and marginal stability curves 
of circular localized states in a rescaled param eter s pace, as 
deduced from the equilibriu m co ndition of Eq. |4.34| and the 
linear stability result of Eq. 



5.22 



FIG. 10. Energy of circular domains as a function of radius. 



from Eq. 4.36, for various values of the activator- inhibitor 
coupling p, with D = 0.01, and r = 0.6. 



FIG. 14. Contour plots of the activator it from numerical 
simulations of the reaction-diffusion model in two dimensions. 
Parameters are D = 0.01, r — 0.40, and p — 0.20. In analogy 
with the one-dimensional front-repulsion phenomenon shown 
in Fig. ^, the approaching chemical fronts do not cross, and 
equilibrate to a simple stripe. Rescaled times are in incre- 
ments of At — 6 from upper left to lower right. 



FIG. 15. Pattern formation with multiple domains of acti- 
vator. D = 0.01, r = 0.60, p = 0.30, Ar = 2.36. The domains 
each undergo fingering instabilities, but do not merge. 



FIG. 16. Numerical simulation of contour dynamics show- 
ing destabilization of a circle. Lower panel shows perimeter 
as a function of time, illustrating the asymptotically linear 
growth. Parameters are: r — 0.212, p = 0.796, with Ar — 25. 



FIG. 17. Tip-splitting in the contour dynamics. Initial con- 
dition is a localized state perturbed by a small n = 2 distor- 
tion. Parameters are f = 0.212, p — 0.90, Ar = 10; both the 
n — 2 and n — 3 modes are unstable. 



FIG. 18. Instability of a 4-fold vertex. Contour dynamics 
results, the initial condition being a 4-fold perturbed local- 
ized state with an additional very small n = 2 distortion. 
f = 0.212, p = 1.10, Ar = 5. 



FIG. 11. Mechanism of instability of a straight chemical 
front. The accumulation of inhibitor in the neighborhood of 
segments of the interface that are concave inwards leads to 
the deepening of the those deformations, and vice- versa for 
regions that are convex outwards. 

FIG. 12. Stability diagram for stripes, in the asymptoti- 
cally rescaled parameter space. The coalescence of the sinu- 
ous instability of a stripe and the instability of a single front 
as r ^ 1/2 refiects the diverging equilibrium width of the 
stripe. Note that the sinuous instability lies in the region of 
azimuthally-stable localized states. 



FIG. 19. Growth of a labyrinth from a compact initial con- 
dition. Contour dynamics results proceed in time from up- 
per left to lower right, spaced by Ar = 5. Parameters are 
f — 0.21, p = 1.10, and initial condition is a distorted circle 
of radius R = 6.0. 

FIG. 20. Relaxation of a labyrinthine pattern to a localized 
state. Initial condition is close to the final panel in Fig. |l^, 
with parameters f = 0.21, p = 0.601, and time between figures 
of Ar = 25. 
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